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Abstract 



We introduce new versions of lattice Boltzmann methods (LBM) for in- 
compressible binary mixtures where fluctuations of total density are inhibited. 
As a test for the improved algorithms we consider the problem of phase sepa- 
ration of two-dimensional binary mixtures quenched from a disordered state 
into the coexistence region. We find that the stability properties of LBM 
are greatly improved. The control of density fluctuations and the possibility 
of running longer simulations allow a more precise evaluation of the growth 
exponent. 

KEY WORDS: Lattice Boltzmann methods; binary mixtures; phase separa- 
tion. 
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I. INTRODUCTION 



Lattice Boltzmann methods (LBM) are computational schemes for solving the macro- 
scopic equations of fluid systems [1,2]. In the case of binary mixtures they have been largely 
used to study the dynamics of systems described by the Navier-Stokes and convection- 
diffusion equations [3]. The advantages of LBM compared with standard methods for par- 
tial differential equations are particularly relevant in problems involving complex boundaries 
such as multiphase flow in porous media [4,1,2]. The large number of applications of LBM in 
studies on complex fluids [5-7], polymer mixtures [8], liquid crystals [9], etc. demonstrates 
the versatility of the method. In studies of spinodal decomposition of binary mixtures a 
version of LBM based on a free-energy approach, where the equilibrium distributions can be 
consistently chosen with the thermodynamics of the system, has been largely used [10-12]. 

The stability properties of LBM are always crucial, especially in phase separation studies 
where long lasting simulations are needed to establish the growth properties. Indeed, the 
average size of domains of the two phases generally grows with power law behavior and one 
wants to determine the exponent [13,14]. However, previous studies showed the existence 
of a stability problem: total density fluctuations can suddenly grow in an imcontrolled way 
interrupting simulations and making problematic the determination of the growth exponent 
[12]. 

In this paper we consider new versions of LBM where the stability of the algorithm 
is improved. Wc have checked the stability properties by studying the problem of phase 
separation of two-dimensional binary mixtures. This is a severe test for stability since 
density fluctuations are enhanced by the presence of many interfaces in the system though 
they should be negligible in incompressible fluids. We flnd that, with the new algorithms, 
density fluctuations are largely suppressed and a more reliable determination of the growth 
exponent can be obtained. 

The incompressible Navier-Stokes equation can be obtained from lattice Boltzmann equa- 
tions if density fluctuations are negligible [1]. Compressibility effects, intrinsic to LBM, can 
produce some serious errors in numerical simulations of incompressible fluids. Due to these 
reasons improved versions of LBM have been introduced for the case of a simple fluid [15,16]. 
New schemes are based on the explicit elimination in the equilibrium distributions of higher 
order terms in the Mach number [17] or on the introduction of a feedback mechanism [18]. 
With these improvements it has been shown that analytical solutions known for particular 
flows are better reproduced than with old LBM versions. 

The above methods are here developed and introduced for the flrst time in the case of 
binary mixtures. The algorithms arc presented in Section II. In Section III wc study the 
problem of phase separation comparing the results obtained with different LBM versions. 
The role of boundary conditions on the problem of stability is examined in Section IV. 
Finally, in Section V, we discuss our results and draw our conclusions. 

II. THE MODEL 

Our simulations are based on the free-energy approach developed by Orlandini et al. 
[10] and Swift et al. [11]. In this scheme a suitable free energy is introduced to control the 
equilibrium properties of the system. 
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A. The lattice Boltzmann scheme 



The free-energy functional used in the present study is 



dr 



^nlnn + |(/^2 + ^/+^(Vv?)2 



(1) 



where n is the total density and ip is the density difference between the two fluids. The term 
in n gives rise to a positive background pressure. The terms in (/? in the free-energy density 
/(n, ip) correspond to the usual Ginzburg-Landau expression typically used in coarse-grained 
models of phase separation [14] . The polynomial terms relate to the bulk properties of the 
fluid. While the parameter b is always positive, a allows to distinguish between a disordered 
(a > 0) and a segregated mixture (a < 0) where two pure phases with (p — ±^J—a/b 
coexist. We will consider quenches into the coexistence region with a < and b = —a 
so that the equilibrium values for the order parameter are (/? = ±1. The gradient term is 
related to the surface tension. The other thermodynamic properties of the fluid follow from 
the free energy (1). The chemical potential difference between the two fluids is given by 
A/i — 5T/5(p — a(p + b(p^ — nV^ip. The thermodynamic pressure tensor is given by 



5n 6(p 

1 
3 



nSap + P^T (2) 



where the chemical part is P^J^"^™ = (2 V^^ + ^^'^'^ ~^ ^'^ (^^'^) ~ 2 ^'^^'^ ^^ap + f^da'pOpip. 

The lattice Boltzmann model is defined on the two-dimensional square lattice with the 
following nine velocity vectors: eo = (0,0), = (cos[(i — l)7r/2], sin[(ii — l)7r/2])c, i = 
1,2,3,4, = (cos[(i - 5)7r/2 + 7r/4], sin[(i - 5)7r/2 + 7r/4])v^c, i = 5,6,7,8, where c = 
Ax/ At, and Ax and At arc the lattice constant and the time step, respectively. Two 
sets of distribution functions fi{r,t) and gi{r,t) evolve according to a single relaxation-time 
Boltzmann equation that is discrete in both time and space [19,20]: 

/,(r + eAt, t + At)- Mr, t) = --[Mr, t) - /^(r, t)], (3) 

r 

gi{r + e,At, t + At)- g,{v, t) = --[^^(r, t) - gl\v, t)], (4) 

where r and T^p are independent relaxation parameters. The distribution functions are 
related to the total density n, to the fluid momentum nv and to the density difference tp 
through 



n 



These quantities are locally conserved in any collision process and, therefore, we require 
that the local equilibrium distribution functions /f^(r, t) and gl'^{r,t) fulfill the Eqs. (5). 
The higher momenta of the local equilibrium distribution functions are defined so that the 
continuum equations pertinent to a binary fluid mixture can be obtained [10,11] 
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Yl /reiaCi^ = C^P^^ + nVaVfi , (6) 



Ygr^ia^^PVa, (7) 



aT^iot^m = c^rA/x5„^ + <^w„u^ . (8) 

where F is a coefficient related to tlie mobility of the fluid. The local equilibrium distribution 
functions can be expressed as an expansion at the second order in the velocity v [10,11]. 

f-'^ = Ai + BjVaeia + Civ"^ + DjVaVfseic^eip + Gi^a^ei^eis z = 1, 2, 3, 4 (9) 

f-"^ = An + BuVaCia + CiiV^ + DjjVaVpeiaCip + Gii^apeiaCip i = 5, 6, 7, 8 

and similarly for the g^'^, i = 0, 8. The coefficients in the previous expansions are explicitly 
written for convenience to compare the original scheme with the alternative ones introduced 
in the following subsections. A suitable choice is [21] 



TDth 1 pthx 

^I,ocp- 1^ , Lr//,a/3-^^ ii4j 

The expansion coefficients for the g'l'^ can be obtained from the previous ones with the formal 
substitutions n ^ tp and TAfiSais- The quantities P^|^ and Afi, which appear in the 

coefficients of the equilibrium distribution functions, can be calculated from (1). 

It has been shown in Refs. [10,11], using a Chapman-Enskog expansion [22], that the 
above described lattice Boltzmann scheme simulates at second order in At the continuity 
equation 

dtn + dainvc,) = + o{At^), (15) 
and the Navier-Stokes equation 
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+31/9, {dtP^f^ - d,{nv^vpv,) - {vpd^P^^'"' + v^d,P;^;n) + o( At^) (16) 

where the isotropic pressure contribution is 

P = cin, (17) 

Cs — c/\/3 and u = {t — l/2)c^At being the speed of sound and the kinematic viscosity, 
respectively. The first fine of Eq. (16) corresponds to the standard Navier-Stokes equation. 
The second hne contains spurious terms. Nonetheless, these terms, being gradient terms 
with higher order derivatives, are negligible compared to the ones in the first line when 
all hydrodynamic fields vary smoothly on the lattice scale [11,12]. Moreover, a convection- 
diffusion equation is also recovered: 

dtip + dc{(pva) = eV^Ai^ + 

-(r^ - l)AtdJ^d^p + %c'Ptr) + «(^^') (18) 

where the macroscopic mobility is = r(r^ — l/2)c^Ai. As before, the spurious terms in 
the second line of Eq. (18) are also expected to be small. 

If the density fluctuations are negligible, the continuity and Navier-Stokes equations 
(15)-(16) describe the behavior of an incompressible fluid. However, this is not always the 
case in LBM where the compressibility effects can produce errors in numerical simulations. 
In the next two subsections we will modify the described LBM in order to reduce unphysical 
compressibility effects. 



B. The limit of small Mach number 

The flrst scheme that we consider and adapt to our LBM is the one introduced in Ref. [17]. 
Previous studies for incompressible fluids have shown that the density fluctuations, 5n, are 
expected to be of order o(M^) in the limit M ^ 0, M — v/cs being the Mach number 
[15,23]. We enforce this condition in the method. By substituting n — Uq + Sn into the 
equilibrium distribution functions (9), Uq being the spccifled constant value of density, and 
neglecting the terms proportional to 5n{'v/c), and 5n{y/cf'^ which are of the order o{M^) 
or higher, we flnd that the net effect is to replace n with tiq into the expansion coefficients 
(11)-(13) keeping unchanged the rest of the scheme. In this way the terms in the equilibrium 
distribution functions (9) are of the order o(M^) or lower. The condition M <g 1 is satisfied 
in our runs where it is always M < 0.01. 

In this case the Chapman-Enskog expansion gives the continuity equation 

daVa^O + o{At'')+o{M^), (19) 

the incompressible Navier-Stokes equation 

dtvp + vad^vp = -—dpp - —daC^P^f"' + vV^p + o{Ae) + o{M^), (20) 

no 77-0 

and the convection-diffusion equation 

dt^ + do,{^v^) = ev'A/x + o{At^) + o{M^), (21) 
where we have not written again the spurious terms of Eqs. (16)- (18). 
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C. A feedback mechanism 



In Ref. [18] a different algoritlim based on a feedback meclianism is proposed to keep tlie 
density nearly constant in LBM for a single fluid. Here we adapt this scheme to our model 
combining it with the one previously described. We use for the expansion coefficients (10) 
the following form: 



Ao 


= aQfi — 




Aj 


— ajTi + 


^ r>chemc 


An 


4 





Let us note that these expressions coincide with (10) when = 4/9 and aj = 1/9. Here we 
only ask that Oq + 5a/ = 1 to guarantee the total density conservation. We consider at each 
node and at each time step before collision the density n fixing the value of oj as 

ai^o^r- 6(1 - — ) (23) 

where = 1/9 and 6 is a positive constant given by experience. In this paper we generally 
use h — 0.01. The value of h is flexible, but if it is too large the simulation will be unstable. 
The value of oq follows from the relation between oq and a/. These values are then used in 
(22). This choice forces, in the next streaming step, more (less) particles to leave the node 
with n > no {n < Uq), resulting in n — no- 

By using a Chapman-Enskog expansion and neglecting terms of order o(M"^) or higher 
we find the continuity equation (19), the incompressible Navier-Stokes equation (20) where 
now the isotropic pressure contribution is 

p = Sajc^n = cln — 3c^6(no — n) + o{6n^) (24) 

and the convection-diffusion equation (21). For the convenience of description, we label the 
old scheme as A, the scheme described in subsection B as B, the scheme combining scheme 
B and the feedback mechanism as C. 



III. PHASE SEPARATION WITH PERIODIC BOUNDARY CONDITIONS 

After a binary mixture is quenched into the coexistence region and domains of the two 
phases are well established, the typical size of domains generally grows as R{t) t"'. A 
simple scaling analysis of the Navier-Stokes and of the convection-diffusion equations shows 
that three regimes can be found depending on the role played by hydrodynamic degrees of 
freedom. At high viscosity the domain growth is governed by a diffusive mechanism and 
q; = 1/3 [24]. When hydrodynamics becomes relevant, the laws R{t) ~ i or R{t) ~ i^/^ 
are expected depending on whether viscous forces or inertial effects dominate, respectively 
[13,14]. However, in real systems the situation is more complex. The physical mechanism 
responsible for viscous growth is not operating in the two-dimensional case [25] and, indeed. 
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this regime has never been observed in simulations [3]. In this paper, since we address the 
question of stabihty which becomes more critical at low viscosities, we will focus on cases 
corresponding to the inertial regime. 

The average size of domains R{t) can be calculated as the first momentum of the structure 
factor, that is 

^ !ik cCk.t) 

Jdk \k\ C{k,t) 

where 

C{k,t) = {ipik,t)ip{-k,t)) (26) 

and (■) is the average over different histories. 

In Fig. 1 we show the behavior of R(t) for the sets of parameters listed in Table I. Similar 
behaviors have been also observed in many other cases here not reported. We have also not 
reported the worst situations occurring at very low r, when the scheme A gets very soon 
unstable {t ~ 100) and the new schemes are not significantly better. All the cases of Table 
I refer to situations where phase separation has already generated well-formed macroscopic 
domains and a significant evolution of R{t) can be observed. For each set, results obtained 
with the schemes A, B and C are compared. While the scheme A can be used only before 
density fiuctuations start to grow indefinitely, the simulations with the new schemes remain 
stable at all times for most of the cases considered. We have found only few cases (e.g., set 
6 at the lowest considered viscosity in Table I) where the schemes B and C are affected by 
instability; however, the instability with the scheme C occurs when the system is almost 
completely phase separated, later than with the scheme B and much later than with the 
scheme A. 

All the cases shown correspond to the inertial regime and for reference we have plotted 
the line with slope 2/3. The sets 1 and 2 only differ in the value of the mobility F. The main 
difference is that, before reaching the hydrodynamic 2/3 regime, the growth is faster when 
the mobility is higher as it can be observed in Fig. 1. The slowing down of growth at large 
times is due to finite size effects appearing when the size of domains becomes comparable 
with the lattice size. 

In order to understand the effects of density fluctuations, we show in Fig. 2 pictures 
of the system for the parameter set 2 at the time t = 5895 when the scheme A becomes 
unstable. Simulations obtained with the old and the new schemes can be compared using 
the same initial conditions. Differences in the order parameter field can be observed due to 
the different treatment of density fluctuations. When the scheme A is used the density plot 
shows the appearing of a peak which will grow indeflnitely. The presence of the instability 
can be also observed in the middle low region of the concentration fleld as a white stripe. 
In the case with the scheme B fluctuations remain limited around the average value n = 1 
and are larger in correspondence of interfaces. The velocity field also shows an instability 
with the scheme A while with the new scheme regular local fiows can be observed. 

/^From sets 3 and 6 of Fig. 1 we observe that R{t) behaves more regularly when calculated 
with the scheme C. This allows a better determination of the growth exponent. The effects 
of using a bigger lattice size can be seen in Fig. 1 for cases 4 and 5. 
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In recent simulations of three-dimensional systems based on the original algorithm the 
growth exponent was obtained by fitting the results of R(t) with the function f{t) = c{t — d)°' 
in a time interval from the inflexion point to the onset of the instability [12]. This procedure, 
as admitted by the same authors of Ref. [12] and discussed also in Ref. [26], is very sensitive to 
small variations of the initial time considered, so that the evaluation of the growth exponent 
is problematic. We have applied a similar procedure to our results confirming the sensitivity 
to small variations of the initial time considered and finding different values for the exponent 
when the old or the new more stable schemes are applied. For example, for the set 2 with 
the old method, we find a — 0.58 if we start from the inflexion point located t — 3000 
and a — 0.61 if data are fltted from t — 4000 (the run becomes unstable at t = 5985). If we 
use the longer series of data obtained with the scheme C we get a = 0.64 starting from the 
inflexion point or at a later time. For the set 4 the fltting procedure for the scheme A gives 
values of a in the range 0.50 — 0.55 changing the initial time for the fltting procedure, while 
the scheme B gives values in the range a — 0.64 — 0.67. Other cases are also worst because 
sometimes, like for the set 3, a reliable determination of the exponent is not possible with 
the old scheme. 

IV. STABILITY WITH WALLS 

In this Section we discuss how the presence of boundary conditions with walls affects 
the stability of simulations of phase separation. Boundary conditions with walls are present 
in problems where a flow is applied to the system, e. g. shear flow, or in situations with 
geometrically complex boundaries. 

Fluctuations of density are generally larger close to the walls due to the fact that the 
propagation and collision steps involve a reflection on the boundary. Then, if the density is 
larger than the average value in a point close to a wall, the reflection tends to enhance this 
fluctuation. An example of the behavior of density just before the system becomes unstable 
is given in Fig. 3 which shows clearly that the instability appears close to the walls. 

Also for the case with walls we have applied the schemes described in Sec. II for con- 
trolling the density fluctuations. In all the considered cases instabilities in simulations with 
walls occur much before than in simulations with periodic boundary conditions. In Table II 
we have reported the stability results for the original and the new scheme C. The comparison 
shows that in most of the cases the new algorithm can be very useful for running longer 
simulations of phase separation. 

V. CONCLUSION 

In this paper we have considered two-dimensional Lattice Boltzmann methods for incom- 
pressible binary mixtures based on a free-energy functional approach. We have introduced 
mechanisms for controlling the total density fluctuations. These mechanisms were already 
shown useful for reproducing more accurately examples of analytical solutions for simple 
fluid flows. 

In the case of binary mixtures, density fluctuations can induce numerical instabilities that 
make impossible to run long enough simulations as it is needed in problems like spinodal 
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decomposition. We have used the new algorithms to study the phase separation of binary 
mixtures in the inertial regime at low viscosities when the effects of density fluctuations are 
worst. We have found that the stability of the simulations can be largely improved. For most 
of the parameters considered, differently from what happens with the original algorithm, it is 
possible to take under control density fluctuations eliminating the occurring of instabilities. 
This allows to study the properties of phase separating binary mixtures more reliably; in 
particular, a more accurate evaluation of the growth exponent can be given. As a future 
research in this direction, it would be interesting to apply the LBM with improved stability 
to the case of a three-dimensional mixture. 
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TABLES 



set At —a, b k, t T 

"1 1 0.0625 004 2 T 

2 1 0.0625 0.04 2 0.2 

3 1 0.00625 0.004 0.5195 2.5 

4 0.2 1.252 X 10"^ 8 x 10"^ 0.521 40 

5 1 0.00313 0.002 0.542 8 

6 1 0.00125 0.0008 0.5015 10 



TABLE I. The sets of parameter used in Figure 1. 
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At 


—a, b 


K 


T 


r 


t^ 

max 


max 


1 


0.083 


0.053 


4.73 


0.1 


> 140460 


> 146450 


1 


0.0625 


0.04 


2 


0.2 


= 10139 


> 184200 


1 


0.0625 


0.04 


1.1 


0.3 


= 2043 


> 10640 


1 


0.00625 


0.004 


0.575 


4 


= 2983 


> 6000 


1 


0.00625 


0.004 


0.5195 


2.5 


= 148 


> 61000 


0.2 


1.252 X 10"^ 


8 X 10-s 


0.521 


40 


= 41.6 


> 1315 


1 


0.00313 


0.002 


0.542 


8 


= 682 


> 10200 


1 


0.00125 


0.0008 


0.5015 


10 


= 97 


= 2120 



TABLE II. Safe running times for the schemes A and C for 8 sets of parameters. t„iax, where 
the upper label refers to the scheme used, is the longest simulated time. The symbol = means that 
the simulation becomes unstable after that time; the symbol > means that the simulation is still 
stable up to that time. 
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FIGURES 




FIG. 1. R{t) behavior for the sets of parameters Hsted in Table I. The sohd, dashed, and dotted 
hnes correspond to scheme A, B, and C, respectively. Runs with the schemes A and B always start 
from random initial configurations. The scheme C is applied either to simulations from random 
initial conditions (sets 5 and 6) or starting from a configuration slightly before the scheme A fails 
to work. Runs with the scheme C were eventually stopped, though still stable, once the system 
was phase separated. 
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Fig. 2 and Fig3 are in JPG format. Captions of Fig.2 and Fig. 3 are as follows, 



Fig. 2 Configurations of the order parameter ip (a)-(b), the total density n (c)-(d), and 
velocity fields (e)-(f) for the parameter set 2 at the time t — 5895. Figures (a), (c), (e) are for 
the scheme A, while (b), (d), (f) are for the scheme B. In picture (a) it is —2.13 < </? < 1.77 
and in (b) it is -1.06 < < 1.05. 

Fig. 3 Three-dimensional plot of the total density n in a case with walls at time t = 1315 
using the scheme B. The parameters are those of set 4 in Table I. 
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This figure "Fig2.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/cond-mat/0309517vl 



This figure "Fig3.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/cond-mat/0309517vl 



